R.  AM. No. 3795 


D D C 


V, 


' The  Numerical  Calculation  of 
Subcritical  Steady  Potential  Flow 
around  an  Unyawed  Ellipsoi^ 


— I 

By  P.  W.|Duci^ 

Theoretical  Aerodynamics  Unit, 
University  of  Southampton 


Reports  and  Memoranda  No.  3795t 
May  5 


^ Summary 

The  problem  of  a subcritical,  potential,  steady,  three-dimensional  flow  past  an  unyawed  ellipsoid  is 
considered,  using  ellipsoidal  coordinates. 

The  full  equations  of  motion  and  the  exact  body  surface  boundary  condition  are  used  throughout.  Further, 
by  means  of  a simple  transformation  the  entire  flow  field  is  taken  into  the  computation.  A finite  difference 
method,  followed  by  an  iterative  process  is  used  for  the  .solution  of  the  flow  equations. 

Mach  number  distributions  are  given  for  a number  of  examples,  for  the  free-stream  flow  aligned  along  either 
the  major  axis,  or  the  second  major  axis  of  the  ellipsoid.  ^ 
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1.  Introduction 


The  problem  of  calculating  the  exact  (in  the  numerical  sense)  solution  for  two-dimensional,  potential, 
subcritical  flows  past  general  shapes  (with  and  without  lift)  has  largely  been  solved,  for  some  years  now, 
primarily  using  the  work  of  Sells'. 

Now,  however,  the  recent  development  of  powerful  computing  machines  has  introduced  the  possibility  of 
calculating  the  flow  around  three-dimensional,  wing-  and  body-like  shapes. 

Here  we  describe  a numerical  scheme  that  calculates  the  subcritical  flow  around  an  unyawed  ellipsoid.  The 
method  is  (numerically)  exact,  the  full,  non-linear  equations  of  motion  being  used  throughout,  the  entire  flow 
field  is  taken  into  the  computation,  and  the  body  surface  boundary  condition  is  satisfied  as  exactly  as  numerical 
differencing  schemes  will  permit. 

Throughout  we  work  in  ellipsoidal  body  coordinates.  This  enables  us  to  satisfy  the  exact  body  surface 
boundary  condition  easily,  and  at  the  same  time  ensures  a refined  mesh  distribution  in  regions  in  which  the 
solution  is  changing  rapidly.  Further,  by  means  of  a simple  transformation  of  one  of  these  body  coordinates,  we 
are  able  to  bring  the  entire  physical  flow  field  into  a finite  working  field.  The  transformation  also  ensures  that 
there  is  a certain  amount  of  bunching  of  mesh  points  near  the  body,  whilst  further  away  the  point  distribution  is 
sparser. 

There  are  essentially  two  unknowns  to  the  problem — the  velocity  potential  and  the  speed  of  sound,  which 
are  connected  by  the  Bernoulli  energy  conservation  equation.  The  velocity  potential  has  a singularity  at 
infinity  in  the  physical  flow  field,  and  this  contribution  is  subtracted  out  of  the  calculation  (being  a known 
quantity)  for  nuni^rical  purposes. 

Section  2 deals  with  the  ellipsoidal  body  coordinate  system,  and  the  equations  of  motion  in  this  coordinate 
system,  together  with  the  special  mathematical  treatment  required  along  the  singular  lines  encountered  in 
these  coordinates.  Section  3 describes  the  numerical  procedures  required  for  computation,  and  the  results  are 
given  in  Section  4.  There  are  three  independent  solutions  to  the  problem  of  the  unyawed  ellipsoid,  one  for  the 
freestream  velocity  aligned  along  each  of  the  three  body  axes.  Here  we  consider  just  two  of  the  three  cases, 
namely  for  the  flow  aligned  along  the  two  longest  axes  of  the  ellipsoid.  The  conclusions  are  given  in  Section  5. 

2.  The  Ellipsoidal  Coordinate  System  and  the  Equations  of  Motion 
2.1.  Ellipsoidal  Coordinates 

In  this  particular  orthogonal,  curvilinear  coordinate  system,  the  cartesian  coordinates  (x,  y,  r),  using  the 
well-known  Jacobian  elliptic  functions  are  given  in  terms  of  the  three  ellipsoidal  coordinates  by 

y=[o-  + (f ')*]*  cn  f ’ cH ^ ’ 
z sn  dn 

where  for  conciseness 

sn  = sn  (^^,  '/a) 
sn  ^ “ sn  (^',  Vl-cr) 

and  similarly  for  the  two  other  elliptic  functions. 

<7  is  a parameter  that  partly  determines  the  shape  of  the  ellipsoid,  and  for  our  purposes 

0<cr<  1. 


In  general 


O^^'^QO 
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where  K and  K are  the  complete  elliptic  integral  and  the  complementary  complete  integral  respectively,  of 
the  first  kind,  viz. 

de 

J()  [1  -cr  sin^  tfl* 


r.J- 

•'0 


[\-a  sin^ 

de 

[1  -(1  -O’)  sin'  ^ 


The  surfaces  = constant  are  an  ellipsoid,  a hyperboloid  of  one  sheet,  and  a hyperboloid  of  two  sheets 

respectively. 

The  equations  defining  these  surfaces  are 

2 2 2 
acn-C^  asn'C' 

— ^=1. 

{l-o-)cn^{^  dn^ 

For  our  purposes  we  set  ^ ‘ = = constant  to  complete  the  definition  of  the  ellipsoid. 

It  is  convenient  to  introduce  the  curvilinear  metrics  of  the  system.  For  a general  orthogonal  system  the  three 
curvilinear  metrics  /4,  (;  = 1,  2,  3)  are  given  by 


Then  for  the  ellipsoidal  coordinate  system 

A,=B2By.B* 

A2  - BlBy 
Ai  = B\Bi 

where 

5*  = [i+©)T*k+(^‘)T‘ 

fli  = [o-  cn  ^ ( 1 - tr) 

B,  = [asnU^^  + (C'f]^. 

For  a more  detailed  examination  of  general  coordinate  systems,  see  Mangier  and  Murray^ 


/ = 1,2,3. 


2.2  Equations  of  Motion 

The  continuity  equation  in  a general  coordinate  system  may  be  written  in  the  form 

where  are  the  covariant  velocity  components,  J is  the  Jacobian  of  the  system,  and  g"  are  the  contravariant 
metric  ten.sors  of  order  two. 

The  speed  of  sound  a,  is  obtained  from  the  Bernoulli  energy  equation,  viz. 

1 2 

7 + -:q  = constant 

y-1  2 
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where  q is  the  fluid  speed,  i.e. 


q^  = g‘'V.V, 

Remembering  that  for  an  o:  ihogonal  system 

(no  summation) 
g,/  = 0 (/5^y). 


Then 


•f  — ^gng22g3i 
— A I A 2^4  3. 


Further  since  the  flow  is  irrotational,  there  exists  a velocity  potential  4>,  where 


V, 


3d> 


/=  1,2,3 


and  equation  (1)  may  be  written  in  the  form 


1')^  a(<7^) 
2a~  dC‘  d{' 


where  A'  = l/v4,  for  orthogonal  systems. 

If  we  now  consider  the  special  case  of  the  ellipsoidal  coordinate  system,  equation  (2)  becomes 


BlB 


1 /.  2<!>,<1>2  ^ ^ 

ilB^V~  BlBlB*'a^r"  B\BlBtB*-a^^  B]B\V  B]bW)- 

2«I>1<1>3  ^ 2<b;<t»3  ^ , 1 / , 'i'3 

B\B\bW^'^  B\B\bW^'^  BiBlV  BiBla')^'^ 
^BlBlB*^('^)  ^ 2BlBlB*-(~^Bi^) 


(2) 


(3) 


where 


and 


(p  = — = (!> 

‘P.  ‘V,. 

. _ 

% •" 


V'  = 


V,  d>. 


(^,)^  (A,)- 


(no  summation). 


We  now  find  it  convenient  to  introduce  the  transformation 


f = 


B*d(' 

dC' 

(l+(^')')*(o-  + (^')V 
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i.e. 

l, ' = tn  C 

where 


lni'  = tn  (f'.Vl-o-) 


and  for  our  purposes 


= sn('lcnC' 


as 


«oo. 


As  well  as  simplifying  slightly  equation  (3),  this  transformation  brings  the  entire  physical  field  into  a finite 
working  field,  and  so  there  is  no  need  to  impose  arbitrary  outer  boundaries  in  order  to  obtain  a finite  range  for 
the  computation.  This  transformation  will  also  ensure  a denser  distribution  of  mesh  points  near  the  body  in  the 
flow  field,  and  fewer  points  as  infinity  is  approached  in  the  physical  field,  when  the  finite  difference 
approximation  is  made. 

After  some  algebra,  the  continuity  equation  (3)  becomes 


1 f.  ‘t’?  24>,ct)2  ^ 1 <t>2  ^ 

WS  B\B\a~r"  B]B\bW 

2<!>2<I>,,  1 / ‘^*3 

BtB\B\a'^''^  B]Bi\  B]BWr^ 


‘I>7 

r 

2B\bW 

<I); 

r B] 

2B\bW 

[b\b\ 

\ 

2B\bW 

vb\b\ 

i)+ 


:)  + 


JL(q2d2\  . 


^2  * —(R^)  + -^—(R^R 


B]Bi 


(4) 


where 


^ d<t> 

<!>/  =-=r 

af 


<t>,/  = 
<!>„- 


a-«t> 

a(f’)’ 

a"<D 

afar' 


The  boundary  condition  on  the  surface  of  the  body  (i'  = C'  = = ^ tn  (I,  is 

<t>,  = 0 

whilst  as  where  <!>»  is  the  freestream  (undisturbed)  potential. 


(5) 
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2.3.  Singular  Lines  of  Transformation 

We  find  that  there  are  singular  lines,  along  which  the  Jacobian  of  the  transformation  vanishes.  This  occurs  if 
any  of  the  B,  =0  (which  implies  that  A,-i,  >l,+i  = 0). 

Since  we  are  dealing  with  a non-zero  ellipsoid,  then  ^ ' 3=  > 0 and  consequently  B,  and  Bj  do  not  vanish  in 

the  flow  field. 

However  B,  = 0 along  the  lines  ^^  = ±K,  ^^  = ±K'  which  correspond  in  the  cartesian  system  to  the 
hyperbolae 

x^z' 

=1,  y = 0. 

1 —O'  a 

Along  these  hyp>erbolae,  A 2 and  A,  vanish  and  special  attention  is  required  in  order  to  establish  the 
equation  of  motion  and  the  fluid  speed  q. 

In  general 

v,v'=g''V,v, 

-ArArAi- 


For  ellipsoidal  coordinates 


^ B\Bl  B\Bl  B\B\ 


We  now  examine  the  behaviour  of  the  above  expression  along  these  singular  lines,  where  Bj  = 0 and  Vj, 
V,  = 0 (this  follows  immediately  from  the  symmetry  of  the  problem). 

Now  if  M,  = <l>,/(B,..iB,+i),  Ui  = ^iKBiBi)  then 

q^  = u]  + ul  + ul, 

where  u]  presents  no  problem. 

In  order  to  calculate  the  sum  of  the  other  two  velocity  components,  we  consider  a second  orthogonal 
coordinate  system  {{',  7/’),  with  origin  at  the  point  (0,  K,  K')  in  the  original  (^’,  (^)  coordinate  system. 

^ ' remains  unchanged  but  77^  and  77  ’ are  inclined  at  an  angle  <1/  to  and  f ^ respectively  as  shown  in  Fig.  1 . 

Then 

r}^  = U^-K)  COS  i/r  - (^ ’ - K')  sin  (/» 

77 ' = — K)  sin  (A  + (^'^  - K')  cos  lA. 


2^  2_  ^ 

“-’"B?B5  BiBi 


cos  sini^f]’  , [<!>„’  cos  sin  (A]^ 

^ bIbI  BtB; 


Applying  L’Hopital’s  rule  twice  to  equation  (7)  we  obtain 


r 2.  21  2[(<D,s7)^-K(t)„7„7)^] 

+ ,„2^  , r„J.77-'t2l  • 

,'_o  (Bi),^„4a-  + (f  ) J 


Alternatively 
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where 


F 


I 

I 


1 


Now 


<t»V=T-T 

dri 


(B])^W  = 


and  because  equations  (8)  and  (9)  must  be  equivalent  (the  flow  speed  must  be  independent  of  the  direction  of 
approach) 


<I>  ’ ’±<I)  2 2 = 0 


Taking  the  positive  sign  we  obtain 


<I>22  + ~ 0. 


(10) 


Taking  the  negative  sign  we  obtain 


*1^22  *5^33  "t"  2 tan  — 0.  (11) 

However  this  latter  equation  implies  that  the  equation  of  motion  along  the  singular  lines  is  dependent  on  the 
direction  of  approach  of  the  lines,  and  so  must  be  neglected  in  favour  of  equation  (10). 

Consequently  in  order  to  calculate  the  flow  speed  (and  hence  Mach  number)  along  the  singular  lines  we  use 

2,  2 (<^>23)^  + (^22)" 

^ cr(l -o-)[cT  + (^')^] 

(<I>23)^  + (^33)" 

cr(l  - a-)l<T  + 

3.  Numerical  Techniques 

In  its  present  form,  equation  (4)  is  unsuitable  for  numerical  treatment.  Firstly  negative  powers  of  Bi  occur, 
and  as  discussed  previously,  Bi=0  along  the  lines  = ±K,  - ±K'.  However  we  have  she  n in  the  previous 

section  that  along  these  singular  lines,  the  continuity  equation  reduces  to  (10).  Thus  for  mesh  points  situated 
along  these  lines,  (10)  must  be  used  in  place  of  (4). 

Secondly  <I)  has  a singularity  at  = AT'  (^ ‘ = oo).  However  this  singular  contribution  to  is  known,  and  may 

be  subtracted  out  of  the  potential,  leaving  a perturbation  potential  <t>  which  is  finite  everywhere,  i.e. 


<I>  = <l>j  + <^. 


This  singularity  arises  from  the  freestream  flow.  Thus  we  can  set 


(12) 

Alternatively  since  we  know  the  exact  incompressible  solution  4>*  (see  Appendix  A),  then  we  can  put 

+ (13) 


I 

4 
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In  both  cases,  the  partial  differential  equation  in  <f>  is  approximated  by  a (central)  finite  difference  equation, 
where  the  derivatives  at  the  point  (^  + /5i,  mS2,  nSj)  are  given  by 

(4>l)l.m.n  = (<^i+  l,m.n  ~ 4>l-l  ,m,n  )+0(S]) 

Zdi 

{^I2)l.m,n  ~ , » « (<^/+I,rti  + l,n  ~ <^/+ 1 ,n - l,n  “ I ,n  + 1 ,n  l.rt - 1 ,n  ) *^ 0(5  | + ^2) 

45,62 

with  similar  expressions  for  the  other  derivatives  of  0.  5,,  52  and  53  are  the  cell  dimensions  in  the  and 
directions  respectively.  Derivatives  of  d>s  may  be  calculated  analytically. 

The  difference  equation,  at  the  point  (^  + /5,,  mS2,  n5,)  may  be  written  in  the  form 

^n4^l,m,n-  1 ^n4^l,m,n  + 1 (14) 

where  a„  contains  the  contributions  in  the  difference  equation  from  the  other  12  neighbouring  mesh  points. 
Details  of  the  a„,  b„,  c„  and  d„  are  given  in  Appendix  B.  Similar  expressions  to  (14)  may  be  written  for  the  I-  or 
m-  (^  or  C^)  directions. 

It  can  be  shown  (Varga^)  that  if  a matrix  is  diagonally  dominant  (the  absolute  value  of  the  diagonal  element 
is  greater  than  the  sum  of  the  absolute  values  of  the  off-diagonal  elements)  then  the  point  Gauss-Seidel 
method  is  convergent.  Further  if  this  method  is  convergent,  then  for  block  tridiagonal  matrices,  such  as  is 
represented  by  a set  of  equations  (14),  with  /,  m fixed,  n varying,  it  is  quicker  to  solve  for  blocks  rather  than 
individual  points  during  the  iterative  process.  For  this  purpose,  a variant  of  Choleski’s  method  was  used  {see, 
for  example,  Hartree”*)  in  which  the  submatrix  is  first  factorised  into  upper  and  lower  matrices  for  inversion. 

Again,  a similar  process  may  be  applied  for  solution  along  lines  of  /,  n fixed  and  m varying,  or  m,  n fixed  and  / 
varying. 

In  order  to  start  the  iterative  process,  the  initial  guess  used  was  the  exact,  incompressible  perturbation 
solution  d>*,  i.e.  if 


then 


whilst  if 


then 


<1>,  =<l>ao  + <t>* 


<(,  = 0. 


For  the  purpose  of  the  computation,  symmetry  was  invoked  fully  in  order  to  reduce  demands  in  both 
computer  storage  and  time. 

Since  the  velocity  potential  was  the  flow  field  quantity  stored  throughout  the  computation,  the  convergence 
criterion  adopted  was  that  the  iteration  process  was  stopped  when  the  maximum  change  in  potential 
encountered  anywhere  within  the  flow  field  was  less  than  1 x 10““'  during  one  iteration.  This  corresponds  to  a 
change  in  local  Mach  number  of  about  1 x 10 

In  order  to  transform  the  final  iterated  solution  into  cartesian  coordinates,  two  linear  interpolation  routines 
were  used.  One  gave  the  Mach  number  distribution  on  the  surface  of  the  ellipsoid,  whilst  the  other  produced 
lines  along  which  the  Mach  number  remained  constant. 

A number  of  runs  was  carried  out  using  different  mesh  sizes  in  order  to  test  the  effect  of  mesh  size  on  two  of 
the  calculations  and  the  results  are  shown  in  Tables  I and  2.  In  the  first  example,  for  which  the  freestream  is 
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aligned  along  the  major  axis  of  the  body  (Table  1),  the  solution  appears  fairly  insensitive  to  changes  of  cell 
dimensions  in  any  of  the  three  body  coordinates.  In  the  other  example  (Table  2)  for  which  the  freestream  is 
aligned  with  the  second  major  axis,  the  calculation  appears  more  sensitive  to  changes  in  mesh  size,  particularly 
in  the  region  of  the  stagnation  point.  The  calculation  in  general  appears  especially  sensitive  to  changes  in  the 
number  of  points  taken  in  the  coordinate  perpendicular  to  the  body  surface.  It  should  be  noted  that  the 
discrepancy  in  the  region  of  the  stagnation  point  in  the  cases  for  which  the  number  of  points  in  the  direction 
was  decreased,  is  largely  due  to  the  linear  interpolation  routine,  because  of  the  rapid  rise  in  Mach  number  in 
this  direction  in  this  region.  Examination  of  the  uninterpolated  results  shows  closer  agreement.  The  results 
seem  to  indicate  convergence  of  local  Mach  number  with  decreasing  mesh  size. 

Generally  the  numerical  system  was  well  behaved,  and  in  order  to  increase  the  rate  of  convergence, 
over-relaxation  was  applied  to  the  velocity  potential.  If  denotes  the  value  of  <t>i.m.n  found  by  block 

inversion  during  the  Nth  iteration  cycle,  and  the  previous  value,  then  the  new  value  is 

<Pl.rn.n- 


where  w is  the  relaxation  parameter. 

Usually  relaxation  parameters  of  1 • 8 to  1 -9  gave  adequate  stability  and  rate  of  convergence.  Solving  along 
lines  of  varying  f generally  gave  the  br st  stability  and  rate  of  convergence.  For  a typical  calculation,  solving 
along  varying  lines,  20-60  iterations  were  required  to  satisfy  the  convergence  criterion.  Solutions  obtained 
along  different  coordinate  directions  were  numerically  identical. 

Further  setting  c}),  instead  of  <l>s  =d>jc,  had  little  effect  on  the  numerical  system.  Local  Mach 

numbers  differed  by  no  more  than  1 x 10"“. 


4.  Results 

Tables  3-8  and  Figs.  3-6  and  8 and  9 show  the  Mach  number  distribution  on  the  surface  of  a number  of 
ellipsoids  at  various  Mach  numbers.  Figs.  2 and  7 show  the  lines  on  the  body  surface  along  which  the  Mach 
number  is  constant. 

For  the  cases  in  which  = x,  the  constant  speed  of  flow  along  the  second  major  axis  of  the  ellipsoid,  present 
in  the  incompressible  solution,  is  still  present,  to  within  2 x 10'“  in  Mach  number. 

For  the  examples  for  which  d>a,  = y,  the  Mach  number  is  approximately  constant  along  the  major  axes, 
although  it  appears  to  increase  as  the  tips  are  approached,  the  variation  in  the  example  in  Table  4 being  about 
0-6  per  cent.  Further  this  effect  was  seen  to  increase  with  an  increase  in  aspect  ratio  and  also  with  an  increase  in 
freestream  Mach  number. 

Fig.  10  shows  the  variation  of  normalised  flow  speed,  at  the  point  on  the  surface  given  in  cartesian 
coordinates  by  (0,0,  with  the  reciprocal  of  aspect  ratio  (aspect  ratio  = (4/Tr)[H-(^,',)‘]*[o'  + (fo)’]"*).  The 
ellipsoids  all  have  thickness/chord  ratios  of  10  per  cent.  For  comparison  the  exact  incompressible  result  is 
shown,  as  are  three  theoretical  values  obtained  for  a 10  per  cent,  two-dimensional,  elliptic  section.  One  value 
was  obtained  by  Rasmussen  and  Heys''  using  a variational  method,  a second  by  Clapworthy*  solving  for 
potential,  and  the  third  by  the  author  using  Sells'  method' . The  results  of  the  latter  two  calculations  are  almost 
identical,  giving  a value  of  q/qx>  = 1 -200.  As  expected  the  present  three-dimensional  calculations  approach  the 
two-dimensional  values  as  the  aspect  ratio  is  increased,  to  within  the  numerical  accuracy  of  the  four  sets  of 
calculations. 


5.  Conclusions 

We  have  shown  that  a finite  difference  approximation  to  the  full  equation  of  motion,  expressed  in  ellipsoidal 
coordinates,  can  produce  a satisfactory  numerical  solution  for  the  three-dimensional  flow  around  an  unyawed 
ellipsoid. 

The  main  advantages  of  using  this  body  coordinate  system  are,  a simplification  of  the  treatment  of  the  body 
surface  boundary  condition,  and  a refined  mesh  distribution  (in  the  physical  flow  field)  where  the  solution  is 
changing  rapidly.  The  only  real  problem,  namely  the  mathematical  treatment  along  the  singular  lines  along 
which  the  Jacobian  of  the  transformation  vanishes,  appears  to  present  no  serious  difficulties. 

The  numerical  scheme  is  generally  well  behaved,  with  usually  just  20  to  60  iterations  required  for  maximum 
changes  in  velocity  potential  of  less  than  1 x lO  * per  iteration,  for  a 20x21  x 20  mesh.  This  corresponds  to 
changes  of  local  Mach  number  of  about  1 x 10  *.  Computations  usually  take  no  more  than  70  seconds  on  a 
CDC  7600. 
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An  extension  of  the  present  work  is  planned  to  ellipsoids  in  mixed  flows.  This  will  require  the  application  of 
retarded  difference  schemes  (Murman  and  Cole^)  in  the  supersonic  regions,  along  rotated  coordinate  axes 
(Jameson**).  Further  the  fore  and  aft  symmetry  invoked  in  the  present  subcritical  case  will  be  lost,  as  a result  of 
the  appearance  of  shock  waves  in  the  flow. 
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b 
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s 

00 
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APPENDIX  A 


Exact  Incompressible  Solution 

The  continuity  equation,  in  ellipsoidal  coordinates,  for  incompressible  flow  reduces  to 

^ a / 1 


g;  a M aa>\  a-O)  , r<t> 


Consider  first  the  case  for  which  the  freestream  velocity  potential  is  given  by 
Then  since  <I>  = x will  be  a solution  of  (A-1)  we  try  for  another  solution,  of  the  form 
remembering jr  =[l+(^’)']* 

Having  substituted  (A-2)  into  (A-1),  we  obtain,  after  some  algebra 

-V  , 


The  solution  for  this  first  order  ordinary  differential  equation  in  /'  may  be  written  in  the  form 


where  the  limits  of  integration  are  chosen  for  convenience,  so  that 

/(^')^0  as  f'^oo 

whilst  Cl  must  be  chosen  such  that  the  surface  boundary  condition  (a<l>/a^ = 0 is  satisfied. 

(A-3)  may  be  integrated  to  give 

C,  , , 

f(c') =- — me,  - E{6,  vW)] 

1 —O’ 

where  F(d,'J\-a),  and  g(g,  Vl-o-)  are  the  incomplete  elliptic  integrals  of  the  first  and  second  kind 
■espectively,  with  modulus  Vl  -a-  and  argument  0 where 

^ = tan“‘ 


rhe  surface  boundary  condition  gives 


( 1 - 0-)fi>/(l-t-(fi)=)((7-t-(^i)‘) 


[( 1 - 0-)  - Cy(\+(d?)(<r  + (Ch?)(F„  - E„)] 


Fa  = F(0a,  -Jl-a)  1 


£■,,  = ^(^0,  vl-O") 
So  = tan"'  X 

SO  ) 
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A similar  analysis  for  «t>»  = y,  assuming  a perturbation  stilution  of  the  form 


4>  = m')y 


yields 


o 


<7(1  -cr)f 


r / — / — (1-0-)^'  1 

- E(S,  V 1 — V 1 -<7)  — , — - = 1 


n/(1+(^')=)(<7  + (^')’)J 


where 


^ ^ <7(i-<7)<r,',v(i+(<ry)(<7+(ri)') 

{<7(  1 - <7)  - [( E„  - <rE,)f,',v/(T+  U!,fUcr  + -liifil-  <7)]} 

where  Eo  and  Eo  are  as  defined  in  (A-4). 


APPENDIX  B 

Details  of  Difference  Equation 

Here  we  assume  throughout  that  we  are  solving  along  lines  of  varying  and  constant  ^ and  but  similar 
expressions  may  be  obtained  for  solutions  in  the  two  other  coordinate  directions.  TTien,  we  can  write  the 
difference  equations  in  the  form 


1 + C„  + b„  1 = <?, 


(B-1) 


subsciipt  h'%  denoting  values  obtained  directly  from  the  solution  of  (B-1),  i.e.  with  no  relaxation  parameter 
applied.  For  the  A’th  iteration,  away  from  the  singular  lines,  we  have 


d d>3 


_ -t(  ^ \ ^ I , ‘1’?  \ 1 I,  <J>3  N) 


4E'?B:,B:;o'>,ar  ABlBWS.dC^ 

d>)  ‘1’:  ^ ,r^2,  , ^3  S 


) 


fliB;V  B]Bla'}  R]B\BW'^'''^  bIbIbW^'^'^  BtBlBla'^^ 


d>;  / Bf  ^ , <t>,  d , , a , \ 


_d_ 

I O 2 ’ 


(b]bI) 

(B\Bl) 
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<t>f  / B'i  d ^ Bi  d , 


‘t>„  S 
B\Bldi^ 


where  iSi,  S2  and  5j  are  the  cell  dimensions  in  the  i"',  ('  and  ^ ’ directions  respectively,  and  where  for  the  .Vth 
iteration 


(*I*/)/.m.n  (‘I’ll )/.m.n  ■*" 


^•'Vl 


2S, 


(0tl)l.m.n  — + 


s\ 


j.lN-1)  ,(NI  j. 

0l*l.m*l.n  0l*l,m.n  0 1 - l.m  0 1 - \.rr  \ ,n 

45,6, 


with  similar  expressions  for  the  other  derivatives  of  the  point  (io  + /5i,  w5,,  nS^). 

0'^'  is  the  value  of  the  (unknown)  potential  obtained  for  the  preset  (Nth)  iteration  cycle,  whilst  0'’^'"  is 
obtained  from  the  previous  cycle. 

Note  we  use  the  latest  value  of  0i,„,„  available.  As  well  as  reducing  computer  store,  this  also  tends  to  increase 
the  rate  of  convergence  of  the  overall  system. 

Along  the  singular  lines,  the  equation  of  continuity  reduces  to  (1 1),  and  so  the  coefficients  become  simply 
(neglecting  any  symmetry  or  antisymmetry) 


‘'■-k 


-<rk) 

^ 0/.rn  l.n 


a„=-( 


55 


)■ 


With  symmetry  fully  invoked,  we  take  the  coordinates  in  the  intervals 


as 

0€/«/„ 

as 

as 

s,=(K'-ii,)n„ 

52  = A/w„ 

8y  = K'/n„. 

Using  the  body  surface  boundary  condition  (5),  using  reflection,  we  obtain 

0 I. m.n  = <^Lm,n  +25|<I>,, 

where  (Cq-Sx,  mS2,  nSj)  is  a point  inside  the  body. 

Symmetry  arguments  give  the  following: 

for  d>x,  = X : 

01.  l.n  = 0l.\.n 
0l,m^*\,n  0l,>n„~l.n 

0l.m.- I = ~0l.rn.\ 

0l.m.n„*  \ 0l,nt,n„-  i • 
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for  <I>oc  = y: 

4>l.m.-  \ — 

The  matrix  equation  (B-1)  is  now  of  block  tridiagonal  form,  viz: 

C()  bu  0 

d\  c,  bi 

0 dz  C2  bi 

0 d„„-l 

dn„  Cn^ 

Note  that  because  of  the  symmetry  of  the  problem  we  can  set 

do  = 0 
=0. 
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TABLE  1 

Variation  of  Local  Mach  Number  Along  y = 0,  with  Mesh  Size  for  <I> 


istribution  6x21x10  10x11x10  10x21x6 
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TABLE  2 

Variation  of  Local  Mach  Number  Alongx  = 0,  with  Mesh  Size  for  «l»x  = y,  Mx  = 0-80,  o-  = 0-04,  ^ = 0-10 
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Fio.  2.  Lines  of  constant  Mach  number. 


Fio.  4.  Mach  number  on  surface  of  ellipsoid. 
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Fig.  6.  Mach  number  distribution  on  surface  of  ellipsoid. 
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Fig.  7.  Lines  of  constant  Mach  number. 


Fig.  8.  Mach  number  distribution  on  surface  of  ellipsoid 


Fig.  10.  Variation  of  normalized  flow  velocity  at  centre  of  upper  (and  lower)  surface  of  ellipsoid  with 

1 /aspect-ratio. 
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